Phase retrieval in phase contrast imaging

ABSTRACT

A method of obtaining in image of the phase change introduced by an object in penetrating radiation incident on the object includes irradiating the object with penetrating radiation having high lateral spatial coherence, and receiving at least a portion of the radiation at a detector after the radiation has merged from the object and thereby obtaining and storing at least two intensity records for the received radiation each including intensity values at predetermined intervals. These values are utilized to derive a grid of values defining an image of the phase change introduced by the object in the penetrating radiation. The intensity records are obtained at a uniform finite distance after the radiation has emerged from the object, and are for respective different energy distributions of the detected radiation. Apparatus is also disclosed.

RELATED APPLICATIONS

This is a continuation of Ser. No. 09/331,861 based upon Ser. No.09/331,861, filed Sep. 13, 1999, now U.S. Pat. No. 6,226,353, which is anational phase of PCT/AU97/00882 filed Dec. 24, 1997.

FIELD OF THE INVENTION

This invention relates generally to the observation of structuralfeatures of objects utilising penetrating radiation such as x-rays. Inparticular, the invention relates to the derivation of images of thephase change introduced by an object in penetrating radiation incidenton the object, from a two-dimensional intensity record of thepenetrating radiation after it has traversed the object. The inventionmay be extended to retrieve separate phase and absorption data from aset of radiographic measurements.

BACKGROUND ART

The present applicant's international patent publications WO 95/05725(PCT/AU94/00480) and WO 96/31098 (PCT/AU96/00178) disclose variousconfigurations and conditions suitable for differential phase-contrastimaging using hard x-rays. Other earlier disclosures of interest are tobe found in Soviet patent 1402871 and in U.S. Pat. No. 5,319,694.Differential phase-contrast imaging shows great promise for viewing theinternal structure of objects for which traditional absorption-contrastradiography is of limited or no value because of very weak absorptioncontrast. This is the case, for example, with soft tissue within thehuman body.

The practical issue of optimally and efficiently deriving thephase-contrast image for an object from the actual record at thedetector is addressed in two related papers by Nugent et al Phys. Rev.Lett. 77, 2961-2964 (1996); J. Opt. Soc. Am. A13, 1670-82 (1996) andreferences therein. In these papers, it has been demonstrated that withmonochromatic plane-wave x-radiation as in the configurations of WO95/05725 and U.S. Pat. No. 5,319,694, the retrieval of phase informationfrom measurements of the propagation of intensity can be based ontreating the propagation of the modified radiation field whosecharacteristics reflect the phase modifying effects of the object. Atwo-dimensional recording of the intensity of the penetrating radiationafter it has traversed the object is the result of variations in thelocal direction of propagation of the radiation arising from variationsin local refractive index, typically an indication of a boundary orrapid variation in electron density within the object or of a thicknessvariation. The aforementioned articles by Nugent et al utilise atreatment of the propagation of a plane monochromatic electromagneticwave based on Maxwell's equations to derive a transport-of-intensityequation and propose solutions of this equation to derive aphase-contrast image from the intensity record. These suggestedsolutions to the transport-of-intensity equation involve expanding thephase in orthogonal functions. The kind of function chosen depends onthe shape of the sample, and thus Zernike polynomials are adequate for acircular shape whilst a Fourier expansion is most suitable for asquare-shaped sample.

The aforementioned international patent publication WO 96/31098discloses an in-line phase-contrast imaging configuration utilising asubstantially point source and a two-dimensional x-ray imaging detectorspaced from the object. It is demonstrated in the application that, incontrast to previous phase-contrast imaging configurations, a pointsource may be utilised, and moreover that the source may be broadlypolychromatic provided its radiation has high lateral spatial coherence,which in practical terms indicates a maximum source diameter (s)dependent upon the source to object distance (R₁). The larger thesource-object distance or the smaller the source size, the greater thelateral spatial coherence (see Wilkins et al Nature 384 335-8 (1996). Aconsequence of these disclosures in WO 96/31098 is that the approachproposed is more closely related to traditional methods used forabsorption-contrast radiography and should be easier to implement thanearlier proposals. This method of phase-contrast imaging is especiallyadvantageous in the hard x-ray region where the lack of suitable lenselements make other techniques conventionally used in visible light andsoft x-ray microscopy unsuitable.

It is an object of the present invention, at least in one or moreembodiments, to provide a method of obtaining a phase-contrast imagefrom a two-dimensional intensity record where the penetrating radiationsubstantially emanates from a point-like source. In one or moreembodiments, it is a particular objective to provide a method adaptableto the extraction of phase and absorption-contrast information fromradiographic images recorded with a microfocus source which need not behighly monochromatic.

SUMMARY OF THE INVENTION

In certain aspects, the invention involves the concept of obtaining twoor more intensity records at a common finite distance after theradiation has emerged from the object, for respective different energydistributions of the radiation.

In one or more other aspects, the invention entails an appreciationthat, with certain features, a point-like source configuration alsolends itself to an approach based on the solution of a differentialtransport-of-intensity equation, albeit a different one from that usedby others for the plane-wave case.

The invention is especially useful for separating out and retrievingphase information from a typical intensity record (obtained with amicrofocus radiation source) which has both phase-contrast andabsorption-contrast content.

In a first aspect, the invention provides method of obtaining an imageof the phase change introduced by an object in penetrating radiationincident on the object, including:

irradiating the object with penetrating radiation having high lateralspatial coherence;

receiving at least a portion of said radiation at detector means afterthe radiation has emerged from the object and thereby obtaining andstoring at least two intensity records for the received radiation eachincluding intensity values at predetermined intervals; and

utilising these values to derive a grid of values defining an image ofthe phase change introduced by the object in the penetrating radiation;

wherein said intensity records are obtained at a uniform finite distanceafter the radiation has emerged from the object, and are for respectivedifferent energy distributions of the detected radiation.

In its first aspect, the invention further provides apparatus forobtaining an image of the phase change introduced by an object inpenetrating radiation incident on the object, including:

means to provide a source for irradiating an object with penetratingradiation having high lateral spatial coherence; and

detector means for receiving at least a portion of said radiation afterthe radiation has emerged from the object whereby to generate at leasttwo intensity records for the received radiation each includingintensity values at predetermined intervals;

wherein said detector means is arranged for obtaining said intensityrecords at a uniform finite distance after the radiation has emergedfrom the object, and energy characterising means is provided wherebysaid intensity records are for respective different energy distributionsof the detected radiation.

In one embodiment, the respective different energy distributions areobtained by altering the energy spectrum of the radiation irradiatingthe object. This might be achieved, for example, by modifying the outputof the radiation source, or by pre-filter means. In another embodiment,the respective different energy distributions are obtained by providingfor the detector means to provide intensity as a function of energy in acertain energy band or band(s). For this purpose, a two-dimensionaldetector means may be variably wavelength sensitive, or be preceded by avariable filter shutter means. Further alternatively, the imageintensity may be recorded as a function of photon energy for each pixelover a number of ranges of x-ray energy. Advantageously, for enhancedresolution, multiple intensity records may be obtained for respectivemultiple energy distributions of the radiation.

In the simplest case, each energy distribution may be a particularwavelength or photon energy level.

The aforesaid derivation may include solving one or more differentialtransport-of-intensity equations relating the phase at a plane of theobject to the evolution of the intensity distribution along thedirection of propagation, utilising predetermined uniform boundaryconditions. An alternative derivation includes solving Fourier-opticsequations. Others are of course possible for particular configurationsor circumstances.

Preferably, the intensity values also reflect absorption contrast in theobject, and the method further includes utilising these values to derivea grid of values defining an effective pure absorption-contrast image ofthe object.

With the apparatus of the first aspect of the invention, there ispreferably further included a computer program product having a set ofmachine readable instructions which, when installed in a computer havinga suitable operating system and memory means, configures the computer tobe operable to utilise said values to derive a grid of values definingan image of the phase change introduced by the object in the penetratingradiation.

In a second aspect, the invention provides a method of obtaining animage of the phase change introduced by an object in penetratingradiation incident on the object, from one or more two-dimensionalintensity records of penetrating radiation after it has traversed theobject, the radiation being of high lateral spatial coherence whenincident on the object and the or each record being obtained at a finitedistance after the radiation has emerged from the object incorporatingphase-perturbed components within a surrounding field of radiationeither uniformly phase-perturbed or not phase-perturbed, the methodincluding:

storing intensity values from the or each record, at predeterminedintervals;

utilising these values and any predetermined uniform boundary conditionsto derive a grid of values defining an image of the phase changeintroduced by the object in the penetrating radiation, by solving adifferential transport-of-intensity equation relating the phase at anexit plane at the object to the evolution of the intensity distributionalong the direction of propagation.

The invention also provides, in its second aspect, a method of obtainingan image of the phase change introduced by an object in penetratingradiation incident on the object, including:

irradiating the object with penetrating radiation having high lateralspatial coherence;

receiving at least a portion of the said radiation at a detector at oneor more finite distances after the radiation has emerged from the objectincorporating phase-perturbed components within a surrounding field ofradiation either uniformly phase-perturbed or not phase-perturbed, andthereby obtaining and storing intensity values for the receivedradiation at predetermined intervals; and

utilising these values and any predetermined uniform boundary conditionsto derive a grid of values defining an image of the phase changeintroduced by the object in the penetrating radiation, by solving adifferential transport-of-intensity equation relating the phase at anexit plane at the object to the evolution of the intensity distributionalong the direction of propagation.

Preferably, the same intensity values are further utilised to alsoderive values defining an effective pure absorption-contrast image forthe object.

The invention further provides, in its second aspect, an apparatus forobtaining an image of the phase change introduced by an object inpenetrating radiation incident on the object, from one or moretwo-dimensional intensity records of penetrating radiation after it hastraversed the object, the radiation being of high lateral spatialcoherence when incident on the object and the or each record beingobtained at a finite distance after the radiation has emerged from theobject incorporating phase-perturbed components within a surroundingfield of radiation either uniformly phase-perturbed or notphase-perturbed, the apparatus including:

(a) means for storing intensity values from the or each record, atpre-determined intervals; and

(b) means, preferably including a stored program of machine readableinstructions, for utilising said values and any predetermined uniformboundary conditions to derive a grid of values defining an image of thephase change introduced by the object in the penetrating radiation aphase-contrast image for the object, by solving a differentialtransport-of-intensity equation relating the phase at a plane at theobject, e.g. an exit plane, to the evolution of the intensitydistribution along the direction of propagation.

The invention also provides a set of machine readable instructionswhich, when installed in computer means also having suitable operatingsystem software and memory means, configures the computer means foroperation as apparatus according to the preceding paragraph. Theinvention still further provides a storage medium, e.g. a magnetic disk,a CD-ROM or optical storage disk, or an internet server, in which isstored said set of machine readable instructions.

The invention still further provides, in its second aspect, apparatusfor obtaining an image of the phase change introduced by an object inpenetrating radiation incident on the object, including:

a source for irradiating an object with penetrating radiation havinghigh lateral spatial coherence;

a detector for receiving at least a portion of said radiation a finitedistance after the radiation has emerged from the object incorporatingphase-perturbed components within a surrounding field of radiation notphase-perturbed or uniformly phase-perturbed, whereby to generateintensity values for the received radiation at pre-determined intervals;and

computer means, including a stored program of machine readableinstructions, operable to utilise said values and any predetermineduniform boundary conditions to derive a grid of values defining an imageof the phase change introduced by the object in the penetratingradiation, by solving a differential transport-of-intensity equationrelating the phase at an exit plane at the object to the evolution ofthe intensity distribution along the direction of propagation.

Preferably, the computer means is further operable to utilise the sameintensity values to also derive values defining an effective pureabsorption-contrast image for the object.

In one embodiment of the method of the second aspect of the invention,two of the intensity records are obtained at different distances afterthe radiation has emerged from the object, e.g. in two distinct imageplanes. In another embodiment two or more of the intensity records areobtained at a uniform distance, e.g. in a single image plane, fordifferent energy distributions of the incident radiation on the object.In a particular form of the latter embodiment, one or more of saidintensity records are recorded such that image intensity as a functionof photon energy is recorded for each pixel over a number of ranges ofx-ray energy.

The transport-of-intensity equation relevant to the point source casemay be equation (16), or in an alternative form, equation (18),hereunder, and the solution may be by a perturbation method.Alternatively, the equation may be solved numerically, especially whenthe last two terms of equation (16) have similar magnitudes.

As employed herein, the term “penetrating radiation” includes x-rays andneutrons, although the invention is especially useful for x-rayradiation, and may be substantially monochromatic but is more typicallybroadly polychromatic. An application of particular interest is therange 0.5 keV to 1 MeV, e.g. the hard x-ray range 1 keV to 1 MeV. Thephase perturbation by the object may be thought of as a refractioneffect or may more rigorously be viewed as a Fresnel diffraction effect.For an object of finite thickness within a surrounding medium, e.g. ofdifferent refractive index, the phase perturbation will also bedependent on the thickness of the object in the direction of thelocalised wave vector.

The object may, e.g., be a boundary, typically a boundary exhibiting asharp variation in refractive index with respect to the penetratingradiation. The invention is especially useful where there is weak ornegligible absorption contrast for the radiation between the intensitiesof the radiation passing through respective sides of the boundary, butmay generally also be applied where there is significant absorptioncontrast at the boundary.

The boundary conditions typically do not need to be measured and may,e.g., include uniform Dirichlet, Neumann or periodic boundaryconditions. The boundary conditions are selected to achieve a uniquesolution of the equation for phase, at least up to an arbitrary constantcomponent.

Preferably, the solution further utilises one or more opticalconditions. Such conditions may include e.g. a small wavefront curvaturefor the incident radiation, absence of focal points between the objectand image, and uniform illumination of the object.

The incident penetrating radiation is not restricted to beingmonochromatic. For polychromatic radiation, the equation includes aspectrally weighted term or factor dependent on the square of therespective wavelength components.

BRIEF DESCRIPTION OF THE FIGURES

The invention will now be further described, by way of example only,with reference to the accompanying figures, in which:

FIG. 1 is a diagram of an x-ray optics configuration according to anembodiment of the invention in which two intensity records are made atdifferent detection planes;

FIG. 2 is a diagram of a related co-ordinate system for the mathematicaldiscussion which follows;

FIG. 3 is a diagram of an alternative x-ray optics configuration inwhich two intensity records are made at a common detection plane forrespective different energy distributions of the radiation;

FIGS. 4 to 9 illustrate successive intensity records in a mathematicaltest case of a method utilising one approach to the invention; and

FIGS. 10 to 15 are photomicrographs depicting a test of an alternativeapproach.

PREFERRED EMBODIMENTS

The configuration illustrated in FIG. 1 includes a microfocus source Sof high lateral spatial coherence and an x-ray two-dimensional imagingdetector D, for example film, photostimulable phosphor plates (eg. Fujiimage plates), or a two-dimensional electronic detector such as acharge-coupled device (CCD) array.

The term “lateral spatial coherence” herein refers to the correlation ofthe complex amplitudes of waves between different points transverse tothe direction of propagation of the waves. Lateral spatial coherence issaid to occur when each point on a wavefront has a direction ofpropagation which does not change over time. In practice, high lateralspatial coherence may, for example, be achieved by using a source ofsmall effective size or by observing the beam at a large distance fromthe source. Generally, lateral coherence length d⊥=λR₁/s, where λ is thex-ray wavelength, R₁ is the source to object distance, and s is themaximum source diameter. For example, for 20 keV x-rays and a source toobject distance of 200 mm, a source size of approximately 20 μm diameteror less would typically be appropriate. The smaller the source size thebetter for the purposes of this invention, provided total flux from thesource is sufficient. Lateral spatial coherence may need to be preservedby careful selection of the x-ray window of the source, e.g. such thatit is of highly uniform thickness and homogeneity. The effects ofpartial spatial and temporal coherence on image contrast and resolutionare described in Pogany et al [Rev. Sci. Instrum. 68 2774-82 (1997)].

Regions of refractive index variation transverse to the local directionof propagation, or thickness or density variations in the direction ofpropagation, can lead to a significant change in the local direction ofpropagation of the wave front passing through those regions.

Thus, referring again to FIG. 1, a spherical wavefront W1 emanating fromthe point source S becomes distorted to W2 on passing through the objectO. By recording the intensity of the wavefront at a sufficient distancefrom the sample, intensity variations due to rapid refractive index andthickness or density variations in the sample may be detected and theirlocation recorded in an image. This corresponds to a form ofdifferential phase-contrast imaging. The location of the imagingdetector is chosen such that the spatial resolution of the detector issufficient to resolve the intensity differences arising from the severedistortions of the wavefront and to optimise contrast subject topractical considerations. The values recorded at the detector alsoreflect absorption in the object and so the intensity values recorded atthe detector also contain extractable absorption-contrast information.

For reasons to be explained subsequently, it is noted here that theobject is disposed so that radiation which emerges from the objectincorporates phase-perturbed components within a surrounding field ofradiation not phase-perturbed or uniformly phase-perturbed.

Typically, the sharp gradients in refractive index or thickness will beimaged as sharp losses or rapid variations in intensity at correspondingpoints in the image. This feature of intensity loss or rapid variationat a given point in the image is to a first approximation essentiallyindependent of wavelength and can therefore lead to very sharp contrastvariations in the image even when a polychromatic source is used.

This configuration has the feature that for a circular sourcedistribution, the spatial resolution in the image is the same for alldirections and is essentially determined by the source size. It also hasthe advantage that considerable magnification of the image is possibleand so recording media such as photostimulable phosphor image plates maybe used which have many desirable properties such as wide dynamic rangeand high sensitivity but not high spatial resolution.

Detector D is typically coupled to suitable computing means M such as apersonal computer utilising a 486 CPU, operating at 66 Mhz, and providedwith appropriate memory and applications software.

Computer means M stores in memory a set of intensity values from theintensity record at detector D, at predetermined intervals intwo-dimensions. In this simple case, the detector is a flat planardetector and the intervals are on a uniform square grid. The manner inwhich these values are derived will depend upon the kind of detector:for example, the detector may be of a pixellated structure in which eachpixel is sampled in a two-dimensional scan and the values fed seriallyto the computer memory. Alternatively, the record may be made and storedin the detector and scanned or sampled as required by the computer.

Computer M further includes a control program by which the storedintensity values from the detector record and any predetermined uniformboundary conditions are utilised to derive a further grid of valuesdefining an effective phase-contrast image in a selected plane of theobject, by solving a differential transport-of-intensity equation.

For simplicity of formal development, we now discuss the theoreticalaspects for the case of a monochromatic point source, with a view todefining a simple differential equation, and solution, applicable to aweak phase object, as described by conditions (4)-(6). The case of abroadly polychromatic source will be considered subsequently.

With reference to the co-ordinate system indicated in FIG. 2, consider amonochromatic point source with wavelength λ=2π/k located at pointz=−R₁, R₁>>λ, on the optical axis Z. The source illuminates an objectwhich is located between the source and the plane (x, y, 0) orthogonalto the optical axis at point z=0, the object being bounded on the righthand side by the plane z=0 (FIG. 2). The object is assumed to be thin,i.e. its z-dimension is much smaller than R₁.

Let $\begin{matrix}{{u_{0}\left( {x,y} \right)} = {\frac{\exp \left\{ {\quad k\sqrt{\left. {x^{2} + y^{2} + R_{1}^{2}} \right\}}} \right.}{\sqrt{1 + {\left( {x^{2} + y^{2}} \right)/R_{1}^{2}}}}{a\left( {x,y} \right)}\exp \left\{ {\quad k\quad {\psi \left( {x,y} \right)}} \right\}}} & (1)\end{matrix}$

be the scalar complex amplitude at the object plane z=0, whereφ(x,y)=kψ(x,y) is the phase “aberration” appearing due to thetransmission of the wave through the object and I(x,y)=a²(x,y) is thecorresponding intensity.

The propagation of the amplitude u(x,y,0)=u₀(x,y) into the source-freehalf-space z>0 is governed by the Helmholtz equation

(∇² +k ²)u=0.   (2)

A solution to the Helmholtz equation, such that u(x,y,0)=u₀(x,y), isgiven by the first Rayleigh-Sommerfeld integral. This integral can bewritten for z>>λ as $\begin{matrix}{{{u\left( {x^{\prime},y^{\prime},z} \right)} \cong {{- \frac{\quad k}{2\pi}}{\int{\int{\frac{\exp \left\{ {\quad k\quad r} \right\}}{r}\frac{z}{r}{u_{0}\left( {x,y} \right)}{x}{y}}}}}},} & (3)\end{matrix}$

where r={square root over ((x′−x)²+(y′−y)²+z²)}.

In order to simply the analysis, we assume that our object is smallcompared to its distance from the source and is weakly scattering, i.e.

ψ(x,y)≡0 and a(x,y)≡1, when (x ² +y ²)≧d ² , d<<R ₁;   (4)

|∂_(x) ^(m)∂_(y) ^(m)ψ(x,y)|<<(R ¹)^(1−(n+m)) , m ² +n ²>0,   (5)

|a(x,y)|≧C>0 and |∂_(x) ^(m)∂_(y) ^(n) a(x,y)|<<Ck|∂ _(x) ^(m)∂_(y)^(n)ψ(x,y)|, m ² +n ²>0   (6)

where R′=R₁R₂/(R₁+R₂) and z=R₂ is the position of the “image” plane,R₂>>λ.

Condition (4) requires the wavefront incident on the object to have asmall curvature. Condition (5) also has a simple physical meaning: itensures, in particular, that there are no focal points between theobject plane z=0 and the image plane z=R₂ (as |∇²ψ|<<1/R′). In theabsence of condition (5) phase distribution in the image plane mayrelate to intensity distribution in the object plane in a more complexmanner due to possible focal phase shifts. Condition (6) essentiallyrequires the incident beam to illuminate the object uniformly and thevariations of the imaginary part of the refractive index inside theobject to be much weaker than the variations of its real part at a givenx-ray wavelength λ. The last requirement applies both for “phase”objects, and also for mixed phase/absorption objects.

Conditions (4 to 6) guarantee that the diffraction angles are small.Therefore, under these conditions the image can be calculated using theFresnel integral: $\begin{matrix}{{{u\left( {x^{\prime},y^{\prime},R_{2}} \right)} = {{- \frac{\quad k}{2\pi}}\frac{\exp \left\{ {\quad k\quad R} \right\}}{R_{2}}{\int{\int{\exp \left\{ {\quad k\quad {S\left( {x,y} \right)}} \right\} {a\left( {x,y} \right)}{x}{y}}}}}},} & (7)\end{matrix}$

where R=R₁+R₂ and $\begin{matrix}{{S\left( {x,y} \right)} = {\frac{\left( {x^{\prime} - x} \right)^{2} + \left( {y^{\prime} - y} \right)^{2}}{2R_{2}} + \frac{x^{2} + y^{2}}{2R_{1}} + {{\psi \left( {x,y} \right)}.}}} & (8)\end{matrix}$

Applying the stationary phase formula to integral (7) we obtain:$\begin{matrix}{{{u\left( {x^{\prime},y^{\prime},R_{2}} \right)} = {{\sum{\frac{\exp \left\{ {\left\lbrack {{kR} + {{kS}\left( {x_{s},y_{s}} \right)} + {\left( {\pi/2} \right){sgn}\quad S_{s}^{''}}} \right\rbrack} \right\}}{R_{2}\sqrt{{\det \quad S_{s}^{''}}\quad}}{a\left( {x_{s},y_{s}} \right)}}} + {O\left( k^{- 1} \right)}}},} & (9)\end{matrix}$

where the sum is over all stationary points (x_(s),y_(s)) correspondingto image point (x′, y′), sgn S″_(s) is the number of negativeeigenvalues of the matrix S″ of second-order partial derivatives of thephase function S evaluated at point (x_(s),y_(s)) and det S″_(s) is itsdeterminant, |O(k⁻¹)|≦const/k. Conditions (6) ensure that the residualterm O(k⁻¹) is much smaller than the leading terms in (9).

Calculating partial derivatives of the phase function (8) and equatingthem to zero, we obtain the following equations for the stationarypoints (x_(s),y_(s)):

x′=Mx _(s) +R ₂ψ′_(x)(x _(s) ,y _(s))

y′=My _(s) +R ₂ψ′_(y)(x _(s) ,y _(s))   (10)

where M=(R₁+R₂)/R₁

and where we used notation of the type ψ_(q)′ for partial derivatives:Ψ_(q)′=∂Ψ/∂q. Equations (10) describe ray trajectories originating onthe object plane at points (x_(s),y_(s)) and ending up at a given point(x′, y′, R₂) of the image; constant M is the coefficient of imagemagnification. Condition (5) ensures that equations (10) have only onesolution, (x_(s), y_(s)), i.e. there is only one ray connecting eachpoint of the image with the corresponding unique point of the object.Consequently, the summation sign in (10) can be omitted.

It is straightforward to calculate the matrix of second-order partialderivatives of the phase function: $\begin{matrix}{{S^{''}\left( {x,y} \right)} = {\begin{pmatrix}{\left( {1/R^{\prime}} \right) + \psi_{xx}^{''}} & \psi_{xy}^{''} \\\psi_{xy}^{''} & {\left( {1/R^{\prime}} \right) + \psi_{yy}^{''}}\end{pmatrix}.}} & (11)\end{matrix}$

This matrix describes the variation of the curvature of the wavefront.In view of condition (5) this matrix is non-degenerate, both itseigenvalues are positive (hence, sgn S″=0) and its determinant can beapproximated by

det S″=(1+R′∇ ²ψ)/R′², R′|∇²ψ|<<1.  (12)

Now we can simplify the stationary phase formula (9): $\begin{matrix}{{u\left( {x^{\prime},y^{\prime},R_{2}} \right)} \cong {\frac{a\left( {x_{s},y_{s}} \right)}{M\sqrt{1 + {R^{\prime}{\nabla^{2}\psi}}}}\exp {\left\{ {\quad {S\left( {x_{s},y_{s}} \right)}} \right\}.}}} & (13)\end{matrix}$

The corresponding intensity distribution on the image plane is

I(x′,y′,R ₂)≅M ² I(x _(s) ,y _(s))[1−R′∇ ²ψ(x _(s) ,y _(s))].  (14)

Let us now use the stationary-point equations (10) to further simplyequation (14): $\begin{matrix}\begin{matrix}{{I\left( {x_{s},y_{s}} \right)} = \quad {I\left( {{{M^{- 1}x^{\prime}} - {M^{- 1}R_{2}\psi_{x}^{\prime}}},{{M^{- 1}y^{\prime}} - {M^{- 1}R_{2}\psi_{y}^{\prime}}}} \right)}} \\{\cong \quad {{I\left( {{M^{- 1}x^{\prime}},{M^{- 1}y^{\prime}}} \right)} - {{R^{\prime}\left( {{\nabla I} \cdot {\nabla\psi}} \right)}{\left( {{M^{- 1}x^{\prime}},{M^{- 1}y^{\prime}}} \right).}}}}\end{matrix} & (15)\end{matrix}$

Substituting (15) into (14) and denoting x=M⁻¹x′, y=M⁻¹y′, we obtain

M ² I(Mx, My, R ₂)≅I(x,y)[1−R′∇²ψ(x,y)−R′∇log I(x,y)·∇ψ(x,y)].  (16)

This equation describes the image of a weak phase object (as describedby conditions (4)-(6)) as a function of the source-to-object andobject-to-image distances. It shows that in the case of negligibleabsorption the image contrast is proportional to “focal” distance R′ andto the Laplacian of the phase shifts introduced into the incident beamon passing through the object. Note that formula (16) in the case of∇ψ=0 and ∇²ψ=0 trivially gives the exact value of intensity in theabsence of any objects between the source and the image plane.

For x-ray wavelengths away from absorption edges of the material of theobject, the linear attenuation coefficient μ is known to beapproximately proportional to λ³ and $\begin{matrix}{{{\psi \left( {x,y} \right)} = {\frac{\lambda^{2}r_{e}}{2\pi}{\int_{- \infty}^{0}{{\rho \left( {x,y,z} \right)}\quad {z}}}}},} & (17)\end{matrix}$

where r_(e) is the classical electron radius and ρ is the electrondensity of the object. These facts and equation (16) imply that in theaccepted approximation the image contrast has a simple quadraticdependence on the wavelength. Therefore, polychromaticity of the sourceis not necessarily an obstacle to the application of this method and canbe explicitly taken into account: the factor λ² is then replaced by aspectrally weighted sum dependent on the square of the respectivewavelength components.

Equation (16) can be used for explicit retrieval of the phase fromintensity measurements in two image planes, indicated at D₁ and D₂ inFIG. 1. To see how this may be achieved, it is convenient to rewrite theequation in the following form:

−∇·[I(x,y)∇ψ(x,y)]=F(x,y),  (18)

where

F(x,y)=[M ² I(Mx, My, R ₂)−I(x,y)]/R′.

Equation (18) relates the transverse derivatives of the phase on theobject plane to the intensity distributions on the object and imageplanes.

Equation (18) is similar to the plane-wave transport-of-intensityequation (TIE) which can be formally obtained from (18) when R₂→0. Thereare, however, some important differences. First, unlike the TIE,equation (18) does not assume the distance between the two planes D₁,D₂where intensity is measured, to be infinitesimally small. Second,equation (18) takes image magnification explicitly into account, and,therefore, is more appropriate for phase imaging with a point source. Itcan be readily seen that an attempt to directly apply the plane-wave TIEto imaging with a spherical incident wave leads to a different functionon the right-hand side of equation (18), namely the distance R′ in thedenominator of F(x,y) is replaced by R₂. The resulting error is due tothe extra wavefront curvature produced by the incident spherical waveand is negligible only when R₂<<R₁. Finally, our method of derivation ofequation (18), which was based on the application of the stationaryphase formula to the Fresnel integral, allowed us to formulate thevalidity conditions (4)-(6) in terms of the amplitude distributionu₀(x,y) on the object plane. These are effectively conditions on theoptical properties of the source, the sample and the geometricalparameters of the imaging layout. In contrast to this, the conventionalmethod for deriving the plane-wave TIE from the Helmholtz equationinvolves requirements on the diffracted amplitude u(x,y,z), z>0, whichare difficult to verify in practice.

In order to find a unique solution to (18) one needs to define someboundary conditions for it. Such boundary conditions can be readilyobtained using condition (4). If Ω is a simply connected domain in theplane (x,y,0) containing the circle of radius d with the centre at theorigin of coordinates, then on the boundary Γ of Ω we have

ψ(x,y)|Γ=0.   (19)

In fact, one can define not only the Dirichlet boundary conditions (19),but any uniform boundary conditions on Γ, e.g. the Neumann or theperiodic ones. This situation is very different from that encountered inastronomical adaptive optics, where the plane-wave TIE is also used forphase retrieval. There the distorted wavefront is not surrounded by theprimary unperturbed or uniformly perturbed one as in the presentsituation and, therefore, boundary conditions must always be measureddirectly.

In the case where the primary unperturbed or uniformly perturbedwavefront is not measured, a procedure of artificially blurring theimage at the edges of the image field or of iteratively retrieving thephase and recalculating the intensity distribution in the region aroundthe edges of the object to effectively create uniform boundaryconditions might be used.

Equation (18) with boundary conditions (19) has a unique solution forphase, provided I(x,y)≧C²>0 in Ω, where C is a constant defined in (6),while solutions to similar problems with the Neumann and periodicboundary conditions are unique up to an arbitrary additive constant. Therequirement for the intensity to be non-zero is satisfied due toconditions (6).

Solutions to the problem (18)-(19) and similar problems with the Neumannand periodic boundary conditions are stable with respect to errors inthe right-hand side function F(x,y). Note, however, that errors inM²I(Mx, My, R₂) and I(x,y) should preferably be small compared not onlyto these intensities, but to their difference which has to be muchsmaller than the intensities themselves according to equation (16) andcondition (5). Therefore, it is preferred, strongly, that the measuredintensity data be fairly “clean” from noise to optimise the reliabilityof quantitative retrieval of the phase.

In the case of uniform intensity, I(x,y)≡I₀ on the object plane a simplesolution to equation (18) with the periodic boundary conditions in thesquare (−D,D)×(−D,D), D>d, can be given by the formula $\begin{matrix}{{\psi_{mn} = {\left( \frac{D}{\pi} \right)^{2}\frac{{\overset{\sim}{F}}_{mn}}{\left( {m^{2} + n^{2}} \right)}}},{{m^{2} + n^{2}} > 0},} & (20)\end{matrix}$

where ψ_(mn) and {tilde over (F)}_(mn) are the Fourier coefficients ofψ(x,y) and the (experimentally measurable) function {tilde over(F)}(x,y)=[M²I (Mx, My, R₂)−I₀]/(I₀R′), respectively.

In the general case of non-uniform intensity distribution on the objectplane, the last two terms in equation (16) can be comparable inmagnitude. However, in many real situations, where the illuminating wavehas uniform intensity and absorption in the object is moderate, but notnegligible, the last term in equation (16) may be smaller than the rest.In such cases equation (18) can be solved by perturbation methods. Inparticular, it is easy to show that an equation similar to (20), butwith the constant intensity I₀ in the numerator of the function {tildeover (F)} (x,y) replaced by the measured variable intensity distributionI(x,y), can be considered as a first order approximation to the exactsolution of equation (18).

Finally, when the last term in equation (16) has a similar magnitude tothe preceding one, system (18)-(19) can be solved numerically using oneof the well known methods for the Laplace equation on the plane.

It is believed that the described derivation and radiographic method andapparatus provide a practical technique for deriving phase-contrastinformation contained in one or more detected intensity distributions intwo dimensions. The configuration shown in FIG. 1 is straight forwardand akin to configurations traditionally used for absorption-contrastradiography. The disclosed method for deriving the phase-contrast imageis particularly suited to planar detectors of a pixellated or gridconstruction or which are read in an incremental manner, and all themore suitable to the combination of such detectors and digitalcomputers.

Furthermore, by processing the detector values for absorption-contrastby conventional methods, one can simultaneously derive bothphase-contrast and absorption-contrast images from, say, radiographicrecordings in two or more two-dimensional planes D₁,D₂.

Recording of two image intensity distributions at different distancesD₁,D₂ in order to carry out phase retrieval may not always be convenientor practical. For example, the measurement of intensity in two planes atdifferent positions along the optic axis usually requires precisionmovements of the detector system D (indicated by arrow d in FIG. 1)which may be undesirable in real-life applications. A second method ofphase retrieval in the present approach may be achieved by measuring twoor more intensity distributions in an image plane at fixed distance fromthe object but with two or more distinct distributions of incidentenergy, e.g. different wavelengths in a simple case or even moreadvantageously by using a 2-dimensional energy-sensitive detector suchas those based on Cd-Mn-Te. With reference again to FIG. 2, consider aparaxial wavefield in the free half-space z≧0 in an arbitrary state oftemporal and spatial coherence with complex amplitude $\begin{matrix}{{U\left( {r,t} \right)} = {\int_{- \infty}^{\infty}{{U\left( {r,\upsilon} \right)}{\exp \left( {{2}\quad \pi \quad \upsilon \quad t} \right)}\quad {\upsilon}}}} & (21)\end{matrix}$

where v=ck/(2π) is the frequency, k is the wavenumber and c is the speedof light in vacuum. Let S(r,v) be the (power) spectral density, then

S(r _(⊥) ,R,v)−S(r _(⊥) ,v)=−R∇ _(⊥)•[S(r _(⊥) ,v)∇_(⊥)ψ(r _(⊥),v)],  (22)

where kψ(r_(⊥),v)=argU(r_(⊥),v) and we omit z=0 from the list ofarguments of all functions. In the monochromatic case we haveU(r,v)=U(r)δ(v−v₀), S(r,v)=I(r)δ(v−v₀) and the phase kψ(r,v₀) coincideswith kψ(r_(⊥))=argU(r_(⊥)).

The TIE can be obtained by integrating equation (22) over the frequencyrange:

I(r _(⊥) ,R)−I(r _(⊥))=−R∇ _(⊥) •[∫S(r _(⊥) ,v)∇_(⊥)ψ(r _(⊥),v)dv],  (23)

where I(r)=∫S(r,v)dv is the (time-averaged) intensity [14]. As alreadydiscussed, this form of the TIE relates the difference between theintensity distributions in the object and the image planes to the phasederivatives in the object plane.

If spatial variations of the spectral density are much weaker than thecorresponding phase variations, then equation (22) can be simplified:

S(r _(⊥) , R,v)=S(r _(⊥) ,v)[1−R ∇ ² _(⊥)ψ(r _(⊥) ,v)].   (24)

This situation is typical for phase objects, i.e. objects with the realpart of the refractive index varying much stronger than its imaginarypart at the given wavelength. This is the case e.g. for thin biologicalsamples imaged with high-energy x-rays. For such objects conventionalimages based on absorption often have poor contrast, whilephase-contrast images are much more informative [5, 8, 9].

Let us assume that the propagation of x-rays through the sample isparaxial and rectilinear. Then the spectral density and the phase in theobject plane can be expressed as

S(r _(⊥) ,v)=S ^(in)(r _(⊥) ,v)exp{−∫μ(r _(⊥) ,z′,v)dz′},   (25)

ψ(r _(⊥) ,v)=ψ^(in)(r _(⊥) ,v)−∫δ(r ₁₉₅ ,z′,v)dz′,  (26)

where the values with superscript “in” correspond to the wave incidenton the sample, integration is through the sample along the linesparallel to the optic axis, 1−δ is the real part of the refractive indexand μ is the linear absorption coefficient. If all wavelengths presentin the spectrum are sufficiently far away from the absorption edges ofthe sample material, then

μ(r,v)≅σ³μ(r,v ₀),  (27)

δ(r,v)≅σ²δ(r,v ₀).  (28)

where σ=v₀/v and v₀ is some fixed frequency from the spectrum.Substituting equations (25)-(28) into equation (24), we obtain

S(R,v)=S ^(in) (v)exp(−σ³ M ₀)[1−R∇ ² _(⊥ψ) ^(in)(v)+Rσ ²∇²_(⊥)ψ₀],  (29)

where we omitted r_(⊥) from the list of arguments of all functions anddenoted M₀=∫μ(r⊥,z′,v₀)dz′ and ψ₀=∫δ(r_(⊥),z′,v₀)dz′. The last two termsin the square brackets in equation (29) must be much less than one iffor example the paraxial approximation is valid.

If the characteristics of the radiation incident on the sample, i.e.S^(in)(v) and ψ^(in)(v), are known, and the spectral density at a givenfrequency at the image plane, S(R,v), can be measured (e.g. with thehelp of an energy-sensitive detector), then equation (29) can be usedfor phase retrieval from two such measurements at different wavelengths.Indeed, if, for example, S(R,v) and S(R,v₀) have been measured, thenwriting e.g. equation (29) for v and v₀ and excluding exp{−σ³M₀} fromthese equations, it is possible to derive $\begin{matrix}{{{- {\nabla_{\bot}^{2}\psi_{0}}} = \frac{1 - \gamma - {R\left\lbrack {{\nabla_{\bot}^{2}{\psi^{in}(v)}} - {{\gamma\sigma}^{3}{\nabla_{\bot}^{2}{\psi^{in}\left( v_{0} \right)}}}} \right\rbrack}}{R\quad {\sigma^{2}\left( {1 - {\gamma\sigma}} \right)}}},} & (30)\end{matrix}$

where $\begin{matrix}{\gamma = {{\frac{S\left( {R,v} \right)}{S^{in}(v)}\left\lbrack \frac{S^{in}\left( v_{0} \right)}{S\left( {R,v_{0}} \right)} \right\rbrack}^{\sigma^{3}}.}} & (31)\end{matrix}$

The phase ψ₀(r_(⊥)) can be recovered by solving the Poisson equation(30). Note that any uniform boundary conditions can be assigned to ψ₀,if the image is surrounded by the unperturbed primary beam. The solutionto equation (30) with typical uniform boundary conditions is alwaysunique at least up to an arbitrary additive constant. It is obvious fromequations (30) and (31) that the stability of phase recovery isdetermined by the ratio σ=v₀/v of the two frequencies and the ratio ofthe spectral densities in equation (31). In an experiment these valuescan be chosen to be sufficiently different from unity, in which case therecovery of the phase using equation (30) will be stable with respect toerrors in the experimental data. Therefore, this technique may have anadvantage over the earlier-described method of phase recovery fromintensity distribution in the beam cross-sections at two differentpositions along the optic axis. Indeed, in this latter case thez-distance between the cross-sections must be small and the differencebetween the two intensity distributions has to be small too, so thecalculation of the phase Laplacian is a numerically unstable procedureof the “division of zero by zero” type.

As with the earlier embodiments of the invention, an analog of equation(30) suitable for in-line phase imaging with a small polychromaticsource can be derived. The additional curvature appearing in thewavefront due to the quasi-spherical nature of the incident beam can beexplicitly taken into account. Consider the following “quasi-spherical”analog of equation (22):

M ² S(Mr _(⊥) ,R,v)−S(r _(⊥) ,v)=−R′∇ _(⊥)•[S(r _(⊥) ,v)∇_(⊥)ψ(r _(⊥),v)]),   (32)

where M=(R₁+R)/R₁ is the magnification coefficient, R′=(R₁ ⁻¹+R⁻¹)⁻¹ andR₁ is the source-to-object distance. Starting from equation (32) andproceeding exactly as above, it is straightforward to derive the Poissonequation for the phase similar to equation (30), but withψ^(in)(r_(⊥),v) replaced by ψ′^(in)(r_(⊥),v)=ψ^(in)(r_(⊥),v)−r²_(⊥)/(2R₁), R replaced by R′ and γ replaced by $\begin{matrix}{\gamma^{\prime} = {M^{2{({1 - \sigma^{3}})}}{{\frac{S\left( {{Mr}_{\bot},R,v} \right)}{S^{in}\left( {r_{\bot},v} \right)}\left\lbrack \frac{S^{in}\left( {r_{\bot},v_{0}} \right)}{S\left( {{Mr}_{\bot},R,v_{0}} \right)} \right\rbrack}^{\sigma^{3}}.}}} & (33)\end{matrix}$

In a slightly different experimental arrangement without anenergy-sensitive detector it may be possible to perform intensitymeasurements in the image plane with two different monochromaticincident beams with frequencies v and v₀. Then the above equations (30)and (32) still hold with the replacement of all spectral densities bythe corresponding intensities.

It is now proposed to outline a method for phase retrieval frompolychromatic images obtained using only a conventional (notenergy-sensitive) detector. Let us assume that the transversalvariations of absorption are so weak that

exp {−σ³ M ₀(r _(⊥))}≡exp{−σ³{overscore (M)}₀}[1−σ³ M′ ₀(r _(⊥))],  (34)

where {overscore (M)}₀ is the “average” absorption in the sample atv=v₀, while |σ³M′₀(r_(⊥))|<<1 everywhere. The constant {overscore (M)}₀can be determined from the knowledge of total intensity incident on thesample and the total intensity of the image. The weak local variationsM′₀(r_(⊥)) of absorption are assumed unknown.

In the first crude approximation one may choose to neglect the termM′₀(r_(⊥)). Then the following Poisson equation for the phase can beeasily derived from equation (29): $\begin{matrix}{{{- {\nabla_{\bot}^{2}\psi_{0}}} = \frac{{I(R)} - {\int{{{S^{out}(v)}\left\lbrack {1 - {R{\nabla_{\bot}^{2}{\psi^{in}(v)}}}} \right\rbrack}{v}}}}{R\quad {\int{{S^{out}(v)}\sigma^{2}{v}}}}},} & (35)\end{matrix}$

using the notation S^(out)(v)=S^(in)(v)exp{−σ³{overscore (M)}₀}.Therefore, in this case the phase can be found from only onepolychromatic image I(R)≡I(r_(⊥),R)=∫S(r_(⊥),R,v)dv by solving equation(35), providing the properties of the incident radiation, i.e. S^(in)(v)and ψ^(in)(v), are known a priori.

It is now proposed to demonstrate that, if the weak local variationsM′₀(r_(⊥)) of absorption have to be taken into account, the phase can befound from two polychromatic images obtained with the incident beamshaving different spectral composition. Substituting equation (34) intoequation (29) and integrating over the frequency range, we obtain

a−bM′ ₀ −c∇ ² _(⊥)ψ₀ +dM′ ₀∇² _(⊥)ψ₀=0.   (36)

where a=I(R)−∫S^(out)(v)[1−R∇² _(⊥)ψ^(in)(v)]dv, b=−∫S^(out)(v)[1−R∇²_(⊥)ψ^(in)(v)]σ³dv, c=R∫S^(out)(v)σ²dv and d=−∫S^(out)(v)σ⁵dv. Given twoimages, I_(j)(R), j=1,2, measured with incident beams having differentspectral composition (described by the a priori known distributionsS_(j) ^(in)(v) and ψ_(j) ^(in)(v), j=1,2), M′₀ can be eliminated fromthe corresponding pair of equation (36) and obtain $\begin{matrix}{{{- {\nabla_{\bot}^{2}\psi_{0}}} = \frac{{a_{1}b_{2}} - {a_{2}b_{1}}}{{b_{1}c_{2}} - {b_{2}c_{1}} + {a_{2}\left( {d_{1} - {d_{2}{b_{1}/b_{2}}}} \right)}}},} & (37)\end{matrix}$

with all the quantities defined for each of the incident beams exactlyas above, e.g. a_(j)=I_(j)(R)−∫S_(j) ^(out)(v)[1−R∇² _(⊥)ψ_(j)^(in)(v)]dv, j=1,2, etc, similarly, it is possible to solve for M₀ toobtain an effective pure absorption-contrast image.

These expressions have demonstrated that the x-ray phase can berecovered using the Transport of Intensity equation from polychromaticimages obtained at a fixed position along the optic axis. It can be donewith a quasi-plane or a quasi-spherical paraxial incident wave, whosespectral density and the phases of the monochromatic components areknown a priori. The approach would preferably involve a substantiallycomplete characterisation of the source prior to phase retrievalexperiments. The above-described embodiments of the inventiondemonstrate that effective phase- and absorption-contrast images can beobtained, in a variety of ways including any one of the following:

(1) From two images obtained with monochromatic incident beams with twodifferent wavelengths (photon energies).

(2) With a polychromatic incident beam and measurement of the imageintensity as a function of x-ray energy using an energy-sensitivedetector (e.g. based on cadmium manganese telluride), with some binningof the energy values to give sufficient image intensity in at least twoenergy ranges (bins).

(3) With a polychromatic incident beam and measurement of the imageintensity using some energy selectivity in the detector, for example bythe use of energy filters based on the use of absorber foils. Inpractice this might involve using a pair of x-ray films with anappropriate filter foil between them to simultaneously record suitableimage intensity data.

(4) From two polychromatic images obtained with incident beams havingdifferent spectral composition and a conventional, rather thanenergy-sensitive, detector.

In each of these above situations, just one measurement of an image issufficient for phase recovery if the transverse spatial variations ofabsorption are negligible. Otherwise, it is possible to obtain separateeffective phase- and absorption-contrast images from the image intensitydata.

FIG. 3 diagrammatically illustrates a variety of arrangements forderiving plural intensity records for respective different energydistributions. Such measurements can be performed, for example, with anenergy-sensitive detector D_(E), or with a source tunable by adjustablecontrol C to vary its wavelength profile, or with a polychromatic sourceS and interchangeable monochromators or filters either ahead of theobject (F₁) or after the object, e.g. as a shutter F₂ front of thenon-energy sensitive CCD detector D. If absorption in the sample isnegligible, the phase can be found from intensity of only onepolychromatic image, provided the properties of the incident radiationare known a priori.

Other approaches to using two or more different spectral distributions(i.e. alternatives to approaches based on transport-of-intensityequations) may involve deconvolution of the spectral distribution fromthe data or, alternatively, a least squares-type fitting procedure toextract ∇²ψ. Even more sophisticated methods based on Bayesian analysismight also be used. It is proposed to now briefly outline phase andabsorption retrieval by a technique based on a Fourier optics approach.

The treatment of phase retrieval in the spherical-wave case is analogousto the treatment for the plane-wave case and quite straightforward. Thissection is thus confined to the latter.

It may be shown that in the plane-wave case, Cowley's form of theKirchhoff's formula for the x-ray wave function can, with certainsmall-angle approximations, be reduced to: $\begin{matrix}{{\psi \left( {x,y} \right)} \simeq {\frac{}{\lambda \quad R_{2}}{\psi_{0}\left( {x,y} \right)}{\exp \left( {{{- }\quad k\quad x^{2}} + \frac{y^{2}}{2R_{2}}} \right)}*{q\left( {x,y} \right)}}} & (38)\end{matrix}$

and therefore

[ψ(x,y)]≡Ψ(f _(x) ,f _(y))=ψ₀(x,y) exp[iπλR ₂(f_(x) ² +f _(y) ²)]Q(f_(x) ,f _(y)),  (39)

where we adopt the convention that real-space functions are inlower-case and the corresponding frequency-space functions (after takingFourier transforms) are in upper-case, e.g. Q(f_(x)f_(y))=[q(x,y)]. Fora further discussion of the steps involved in this derivation, see J. M.Cowley (1975) Diffraction Physics (North-Holland: Amsterdam.). Relevanttheoretical discussion is also to be found at Pogany et al (1997) Rev.Sci. Inst. 68, 2774-2782.

The transmission function can be written, assuming φ(x,y)t(x,y) andμ(x,y)t(x,y) are sufficiently small, as q(x,y)≅1+iφ′(x,y)−μ′(x,y), whereφ′(x,y)=−φ(x,y)t(x,y) and μ′(x,y)=μ(x,y)t(x,y)/2. The Fourier transformof the transmission function is then given byQ(f_(x),f_(y))≅δ(f_(x),f_(y))+iΦ′(f_(x),f_(y))−M′(f_(x),f_(y)). Equation(39) can then be written as

Ψ(f _(x) ,f _(y))≅ψ₀exp[iX(f _(x) ,f _(y))][δ(f _(x) ,f _(y))+iΦ′(f _(x),f _(y))−M′(f _(x) ,f _(y))]  ( (40)

where X(f_(x),f_(y))=πλR₂(f_(x) ²+f_(y) ²). By expanding (40) and takinginverse Fourier transforms of both sides, one can obtain, to first orderin φ′ and μ′, and omitting the functional dependence:

I≅1−2φ′* ⁻¹[sin(X)]−2μ′* ⁻¹[cos(X)].  (41)

Consequently, in the case of a pure-phase object, one can retrieve theφt-distribution from a single image: $\begin{matrix}{{{\varphi \quad t} = {^{- 1}\left\lbrack \frac{\left( {I - 1} \right)}{2{\sin (x)}} \right\rbrack}},} & (42)\end{matrix}$

and, in the case of a pure-absorption object: $\begin{matrix}{{\mu \quad t} = {- {{^{- 1}\left\lbrack \frac{\left( {I - 1} \right)}{\cos (x)} \right\rbrack}.}}} & (43)\end{matrix}$

In the more usual case of an object for which both phase and absorptioneffects are significant, two images, I₁ and I₂, which have beencollected with different values of R₂ and/or x-ray energy (we arediscussing the monochromatic case here, but in the polychromatic casetwo different voltage settings of the x-ray source, giving two differentspectral distributions, would be appropriate) could be used. Assuming,in the absence of any absorption edges for constituent elements, thatμ₂=(λ₂/λ₁)³μ₁ and that φ₂=(λ₂/λ₁)φ₁, the resulting simultaneousequations can be solved to yield: $\begin{matrix}{{\varphi \quad t} = {- {^{- 1}\left\lbrack \frac{{{\left( {I_{2} - 1} \right)}{\cos \left( x_{1} \right)}} - {\left\lbrack \frac{\lambda_{2}}{\lambda_{1}} \right\rbrack^{3}{\left( {I_{1} - 1} \right)}{\cos \left( x_{2} \right)}}}{2\left\lbrack {{\left\lbrack \frac{\lambda_{2}}{\lambda_{1}} \right\rbrack^{3}{\sin \left( x_{1} \right)}{\cos \left( x_{2} \right)}} - {\left\lbrack \frac{\lambda_{2}}{\lambda_{1}} \right\rbrack {\sin \left( x_{2} \right)}{\cos \left( x_{1} \right)}}} \right\rbrack} \right\rbrack}}} & (44)\end{matrix}$

and $\begin{matrix}{{\mu \quad t} = {{^{- 1}\left\lbrack \frac{{\left\lbrack \frac{\lambda_{2}}{\lambda_{1}} \right\rbrack {\left( {I_{1} - 1} \right)}{\sin \left( x_{2} \right)}} - {{\left( {I_{2} - 1} \right)}{\sin \left( x_{1} \right)}}}{{\left\lbrack \frac{\lambda_{2}}{\lambda_{1}} \right\rbrack^{3}{\sin \left( x_{1} \right)}{\cos \left( x_{2} \right)}} - {\left\lbrack \frac{\lambda_{2}}{\lambda_{1}} \right\rbrack {\sin \left( x_{2} \right)}{\cos \left( x_{1} \right)}}} \right\rbrack}.}} & (45)\end{matrix}$

The advantages of using the multiple energy method of retrieving phaserather than the multiple distance method include the following:

ability to record data by electronic rather than mechanical change ofimaging conditions (for example by rapid switching of tube voltage orenergy-sensitive detectors) leading to the possibility of very rapidphase retrieval; and

ability to use 2-dimensional energy dispersive detectors (e.g. thosebased on CdMnTe) as these become available, in order to record a largenumber of 2-dimensional images for different x-ray energy ranges.

EXAMPLE 1

To test the validity of the method of the invention utilisingtransport-of-intensity equations, a numerical example was prepared withsimulated data, and the successive steps are illustrated in FIGS. 4 to9.

A perfectly monochromatic x-ray point source was assumed with wavelengthλ=0.154 nm (CuK_(a) radiation), and source-to-object and object-to-imagedistances R₁=R₂=0.2 m, so that R′=0.1 m and M=2 (FIG. 2). A pure phaseobject was also assumed, of size 640×640 μm², producing the phasedistribution on the object plane (x,y,0) shown in FIG. 4. The range ofthe phase values was [−0.8667, 0.6706] radians.

The intensity distribution I(x,y,R₂) on the image plane was thencalculated by computing the Fresnel integral (7 above) on a grid with128×128 pixels. The scaled image M²I(x/M,y/M,R₂) of the obtainedintensity distribution I(x,y,R₂) is presented in FIG. 5. For comparisonthe calculated distribution of the phase Laplacian, −∇²φ(x,y,0) wasplotted and is shown in FIG. 6. Obviously, the scaled distribution ofintensity on the image plane and the distribution of the Laplacian ofthe phase on the object plane appear almost identical, which supportsformula (16) in the case of negligible absorption.

At the next stage the computed intensity distribution on the image planeand formula (20) were used to reconstruct the phase. The result of thisreconstruction is shown in FIG. 7. Again, it is in a very good agreementwith the original phase from FIG. 4, with the relative mean-square errorof this reconstruction equal to 4.65%.

The performance of the method was also tested in the case of non-uniformintensity, i.e. when the last term in formula (16) was significant. Anintensity distribution on the object plane was simulated according tothe formula:

I(x,y)=−1.647+3x exp {−(x ² +y ²)/(2*1280²) }.

The values of the intensity had the range of [1, 1.1243]. The intensitydistribution I(x,y,R₂) on the image plane was computed by evaluating theFresnel integrals. The non-uniformity of intensity distribution resultedin the “vignetting” of the corresponding image (FIG. 8), which clearlyreflects the non-uniformity of intensity on the object plane. However,by subtracting the intensity distribution on the object plane from thaton the image plane, dividing the difference by the scaled interplanardistance R′=0.1 m and applying formula (2) we were able to reconstructthe original phase with a root mean-square error of 6.67% (FIG. 9). Thisresult confirms that in the case of moderate non-uniformity of intensitydistribution on the object plane (e.g. weak absorption effects and/orweak non-uniformity of the incident beam), the described simplevariation of formula (20) can be effectively used for quantitativeretrieval of the phase shifts produced by the object.

EXAMPLE 2

FIGS. 10 to 15 show an example of phase and absorption retrieval fromtwo images at different values of R₂. The CSIRO logo has been used asthe absorption (μt) and phase (φt) distributions (FIGS. 10, 11respectively). For the former the values are either 0.0 (black) or 0.1(white), and for the latter (inverted logo) −0.1 (black) or 0.0 (white).The images (intensity distributions) were calculated for the plane-wavecase, 1 Å radiation, 0.2×0.2 mm² (512×512 pixels) area, and R₂=80 cm(FIG. 12) and R₂=160 cm (FIG. 13). These conditions are such that theimages are not immediately interpretable. Applying the retrievalalgorithm (44, 45) discussed, we obtain the results shown (FIGS. 14, 15)which while not perfect, are clearly recognisable.

What is claimed is:
 1. In a method of obtaining an image of the phasechange introduced by an object in penetrating radiation incident on theobject, the method including the steps of irradiating the object withpenetrating radiation having high lateral spatial coherence; andreceiving at least a portion of said radiation at detector means afterthe radiation has emerged from the object, the improvement comprising:obtaining and storing at least two intensity records for the receivedradiation each including intensity values at predetermined spatialintervals; and utilising these values for derivation of a grid of valuesdefining an image of the phase change introduced by the object in thepenetrating radiation; wherein said intensity records are obtained atthe same finite distance after the radiation has emerged from theobject, and are for respective different energy distributions of thedetected radiation.
 2. A method according to claim 1 wherein therespective different energy distributions are obtained by altering theenergy spectrum of the radiation irradiating the object.
 3. A methodaccording to claim 1 wherein the respective different energydistributions are obtained by providing for said detector means toprovide intensity as a function of energy in a certain band or band(s).4. A method according to claim 1 wherein said derivation includessolving one or more differential transport-of-intensity equationsrelating the phase at a plane of the object to the evolution of theintensity distribution along the direction of propagation, utilisingpredetermined uniform boundary conditions.
 5. A method according toclaim 1 wherein said derivation includes solving Fourier-opticsequations.
 6. A method according to claim 1 wherein said intensityvalues also reflect absorption contrast in the object, and the methodfurther includes utilising said values to derive a grid of valuesdefining an effective pure absorption-contrast image of the object.
 7. Amethod according to claim 1 wherein said penetrating radiation comprisesx-ray radiation.
 8. A method according to claim 7 wherein said x-rayradiation is in the range 0.5 KeV to 1 MeV.
 9. A method according toclaim 7 wherein said irradiating radiation is substantiallymonochromatic.
 10. A method according to claim 7 wherein saidirradiating radiation is polychromatic.
 11. A method according to claim7 wherein said irradiating radiation is from a substantially pointsource of full width half-maximum 40 μm or less.
 12. In an apparatus forobtaining an image of the phase change introduced by an object inpenetrating radiation incident on the object, the apparatus includingmeans to provide a source for irradiating an object with penetratingradiation having high lateral spatial coherence, the improvementcomprising: detector means for receiving at least a portion of saidradiation after the radiation has emerged from the object whereby togenerate at least two intensity records for the received radiation eachincluding intensity values at predetermined spatial intervals; whereinsaid detector means is arranged for obtaining said intensity records atthe same finite distance after the radiation emerged from the object,and energy characterising means is provided whereby said intensityrecords are for respective different energy distributions of thedetected radiation.
 13. Apparatus according to claim 12 wherein saidenergy characterising means is arranged to alter the energy spectrum ofthe radiation irradiating the object.
 14. Apparatus according to claim13 wherein said energy characterising means includes means rendering thedetector means able to provide intensity as a function of energy in acertain energy band or band(s).
 15. Apparatus according to claim 12further including a computer program product having a set of machinereadable instructions which, when installed in a computer having asuitable operating system and memory means, configures the computer tobe operable to utilise said intensity values for derivation of a grid ofvalues defining an image of the phase change introduced by the object inthe penetrating radiation.
 16. Apparatus according to claim 15 whereinsaid derivation includes solving one or more differentialtransport-of-intensity equations relating the phase at a plane of theobject to the evolution of the intensity distribution along thedirection of propagation, utilising predetermined uniform boundaryconditions.
 17. Apparatus according to claim 15 wherein said derivationincludes solving Fourier-optics equations.
 18. Apparatus according toclaim 12 further including a source of x-ray radiation as said sourcefor irradiating the object.
 19. Apparatus according to claim 18 whereinsaid x-ray radiation is in the range of 0.5 KeV to 1 MeV.
 20. Apparatusaccording to claim 18 wherein said irradiating radiation issubstantially monochromatic.
 21. Apparatus according to claim 18 whereinsaid irradiating radiation is polychromatic.
 22. Apparatus according toclaim 18 wherein said source is a substantially point source of fullwidth half-maximum 40 μm or less.
 23. In a method of obtaining an imageof the phase change introduced by an objection penetrating radiationincident on the object, from one or more two-dimesional intensityrecords of penetrating radiation after it has traversed the object, theradiation being of high lateral spatial coherence when incident on theobject and the or each record being obtained at a finite distance afterthe radiation has emerged from the object incorporating phase-perturbedcomponents within a surrounding field of radiation either uniformlyphase-perturbed or not phase-perturbed, the method including storingintensity values from the or each record, at predetermined spatialintervals; and utilising these values and any predetermined uniformboundary conditions to derive a grid of values defining an image of thephase change introduced by the object in the penetrating radiation, theimprovement comprising: solving a differential transport-of-intensityequation relating the phase at a plane at the object to the evolution ofthe intensity distribution along the direction of propagation.
 24. Amethod according to claim 23 wherein said intensity values also reflectabsorption contrast in the object, and the method further includesutilising said values to derive a grid of values defining an effectivepure absorption-contrast image of the object.
 25. A method according toclaim 23 wherein said penetrating radiation comprises x-ray radiation.26. A method according to claim 25 wherein said x-ray radiation is inthe range 0.5 KeV to 1 MeV.
 27. A method according to claim 25 whereinsaid irradiating radiation is substantially monochromatic.
 28. A methodaccording to claim 25 wherein said irradiating radiation ispolychromatic.
 29. A method according to claim 28 wherein said equationincludes a spectrally weighted term or factor dependent on the square ofthe respective wavelength components.
 30. A method according to claim 23wherein said boundary conditions include uniform Dirichlet, Neumann orperiodic boundary conditions, and are selected to achieve a uniquesolution of the equation for phase, at least up to an arbitrary constantcomponent.
 31. A method according to claim 30 wherein the solutionfurther utilises one or more optical conditions selected from the groupconsisting of a small wavefront curvature for the incident radiation,absence of focal points between the object and image, and uniformillumination of the object.
 32. In a method of obtaining an image of thephase change introduced by an object in penetrating radiation incidenton the object, the method including the steps of irradiating the objectwith penetrating radiation having high lateral spatial coherence;receiving at least a portion of the said radiation at a detector at oneor more finite distances after the radiation has emerged from the objectincorporating phase-perturbed components within a surrounding field ofradiation either uniformly phase-perturbed or not phase-perturbed, andthereby obtaining and storing intensity values for the receivedradiation at predetermined spatial intervals; and utilising these valuesand any predetermined uniform boundary conditions to derive a grid ofvalues defining an image of the phase introduced by the object in thepenetrating radiation, the improvement in the method comprising: solvinga differential transport-of-intensity equation relating the phase at aplane at the object to the evolution of the intensity distribution alongthe direction of propagation.
 33. A method according to claim 32 whereinsaid intensity values also reflect absorption contrast in the object,and the method further includes utilising said values to derive a gridof values defining an effective pure absorption-contrast image of theobject.
 34. A method according to claim 32 wherein said penetratingradiation comprises x-ray radiation.
 35. A method according to claim 34wherein said x-ray radiation is in the range 0.5 KeV to 1 MeV.
 36. Amethod according to claim 34 wherein said irradiating radiation issubstantially monochromatic.
 37. A method according to claim 34 whereinsaid irradiating radiation is polychromatic.
 38. A method according toclaim 37 wherein said equation includes a spectrally weighted term orfactor dependent on the square of the respective wavelength components.39. A method according to claim 32 wherein said boundary conditionsinclude uniform Dirichlet, Neumann or periodic boundary conditions, andare selected to achieve a unique solution of the equation for phase, atleast up to an arbitrary constant component.
 40. A method according toclaim 39 wherein the solution further utilises one or more opticalconditions selected from the group consisting of a small wavefrontcurvature for the incident radiation, absence of focal points betweenthe object and image, and uniform illumination of the object.
 41. In anapparatus for obtaining an image of the phase change introduced by anobject in penetrating radiation incident on the object, the apparatusincluding a source for irradiating an object with penetrating radiationhaving high lateral spatial coherence; and a detector for receiving atleast a portion of said radiation a finite distance after the radiationhas emerged from the object incorporating phase-perturbed componentswithin a surrounding field of radiation not phase-perturbed or uniformlyphase-perturbed, whereby to generate intensity values for the receivedradiation at pre-determined spatial intervals, the improvementcomprising: computer means, including a stored program of machinereadable instructions, operable to utilise said values and anypredetermined uniform boundary conditions to derive a grid of valuesdefining an image of the phase change introduced by the object in thepenetrating radiation, by solving a differential transport-of-intensityequation relating the phase at a plane at the object to the evolution ofthe intensity distribution along the direction of propagation.